TeaTime4schools: Joint analysis - fungi

07 Beta diversity analysis

roey.angel@bc.cas.cz

Beta diversity analysis

This analysis explores the alpha-diversity ditribution patters in the different samples, based on the DADA2-produced sequences. ### Setting general parameters:

set.seed(1000)
min_lib_size <- 5000
data_path <- "./DADA2_pseudo/"
Ps_file <- "TeaTime4Schools_ITS_filt3.RDS"
Proj_name <- "TeaTime4Schools"

Load phyloseq object

This phyloseq object was created in 05_Taxonomical_analysis.html. The Ps_obj_filt object excludes contaminants or unknown but still includes taxa with low prevalence

readRDS(file = paste0(data_path, Ps_file)) %>%
  subset_samples(., Field != "Unburied") %>% # drop unburied samples
  prune_samples(sample_sums(.) > min_lib_size, .) %>% # remove samples  < min_lib_size
  filter_taxa(., function(x) sum(x) > 0, TRUE) -> # drop taxa with 0 abundance
  Ps_obj_filt

qplot(rowSums(otu_table(Ps_obj_filt)), geom = "histogram") + 
  xlab("Library size")

qplot(log10(rowSums(otu_table(Ps_obj_filt)))) +
  xlab("Logged library size")

Standardize abundances to the median sequencing depth (and convert to proportion)

adonis(
  otu_table(Ps_obj_filt) ~ Lib.size,
  data = get_variable(Ps_obj_filt),
  method = "horn",
  permutations = 9999
)

## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt) ~ Lib.size, data = get_variable(Ps_obj_filt),      permutations = 9999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Lib.size    1      3.57  3.5697  10.438 0.08321  1e-04 ***
## Residuals 115     39.33  0.3420         0.91679           
## Total     116     42.90                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ps_obj_filt %>%
  otu_table(.) %>%
  as(., "matrix") %>%
  rowSums() %>% 
  median() ->
  total
standf = function(x, t = total) round(t * (x / sum(x)))
Ps_obj_filt_median <- transform_sample_counts(Ps_obj_filt, standf)  # Standardize abundances to median sequencing depth

Ps_obj_filt_median_rel <- transform_sample_counts(Ps_obj_filt_median, function(x) x / sum(x)) # convert to relative abundance (just incase it's explicitly needed)

sample_data(Ps_obj_filt_median)$Lib.size <- sample_sums(Ps_obj_filt_median)

qplot(rowSums(otu_table(Ps_obj_filt_median)), geom = "histogram") + 
  xlab("Library size")

PlotLibDist(Ps_obj_filt_median)

PlotReadHist(as(otu_table(Ps_obj_filt_median), "matrix"))

notAllZero <- (rowSums(t(otu_table(Ps_obj_filt_median))) > 0)
meanSdPlot(as(log2(t(otu_table(Ps_obj_filt_median))[notAllZero, ] + 1), "matrix"))

Variance partitioning models and ordinations

Partitioning the data using discrete distance

adonis(
  otu_table(Ps_obj_filt_median) ~ Lib.size,
  data = get_variable(Ps_obj_filt_median),
  method = "horn",
  permutations = 9999
)

## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median) ~ Lib.size, data = get_variable(Ps_obj_filt_median),      permutations = 9999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Lib.size    1     0.122 0.12239 0.32903 0.00285 0.9683
## Residuals 115    42.778 0.37198         0.99715       
## Total     116    42.900                 1.00000

adonis(
  otu_table(Ps_obj_filt_median_rel) ~ Field + Sample.type * Season,
  data = get_variable(Ps_obj_filt_median_rel),
  method = "horn",
  permutations = 999
)

## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_rel) ~ Field +      Sample.type * Season, data = get_variable(Ps_obj_filt_median_rel),      permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                     Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field                2     2.327  1.1636   7.969 0.05425  0.001 ***
## Sample.type          2     9.741  4.8706  33.356 0.22707  0.001 ***
## Season               3     9.463  3.1543  21.602 0.22058  0.001 ***
## Sample.type:Season   6     6.329  1.0548   7.224 0.14752  0.001 ***
## Residuals          103    15.040  0.1460         0.35058           
## Total              116    42.900                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# sample_data(Ps_obj_filt_median)$Type.season <- paste(sample_data(Ps_obj_filt)$Sample.type, sample_data(Ps_obj_filt)$Season) # too many comparisons

# Ps_obj_filt_median_s <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median) # for testing

# Compare sample-type pairs
mod_pairwise1 <- PairwiseAdonis(
  x = otu_table(Ps_obj_filt_median),
  factors = get_variable(Ps_obj_filt_median, "Sample.type"),
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise1)

##                  pairs total.DF   F.Model         R2 p.value p.adjusted sig
## 1 Green tea vs Rooibos       86  7.952102 0.08555053   0.001      0.001  **
## 2    Green tea vs Soil       74 29.422941 0.28726905   0.001      0.001  **
## 3      Rooibos vs Soil       71 18.848044 0.21213797   0.001      0.001  **

(sig_pairs1 <- as.character(mod_pairwise1$pairs[mod_pairwise1$p.adjusted < 0.05]))

## [1] "Green tea vs Rooibos" "Green tea vs Soil"    "Rooibos vs Soil"

# (meaningful_sig_pairs1 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))

# compare seasons
mod_pairwise2 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median),
  sample_data(Ps_obj_filt_median)$Season,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise2)

##              pairs total.DF   F.Model         R2 p.value p.adjusted sig
## 1 Winter vs Spring       55  5.013270 0.08495157   0.001     0.0012   *
## 2 Winter vs Summer       58 22.846552 0.28613073   0.001     0.0012   *
## 3 Winter vs Autumn       61 17.165072 0.22244613   0.001     0.0012   *
## 4 Spring vs Summer       54 14.048243 0.20952440   0.001     0.0012   *
## 5 Spring vs Autumn       57  9.216003 0.14131505   0.001     0.0012   *
## 6 Summer vs Autumn       60  2.448077 0.03983977   0.017     0.0170   .

(sig_pairs2 <- as.character(mod_pairwise2$pairs[mod_pairwise2$p.adjusted < 0.05]))

## [1] "Winter vs Spring" "Winter vs Summer" "Winter vs Autumn" "Spring vs Summer"
## [5] "Spring vs Autumn" "Summer vs Autumn"

# (meaningful_sig_pairs2 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))

Calculate ordinations

Ps_obj_ord1 <- ordinate(Ps_obj_filt_median, "CAP", "horn", formula = Ps_obj_filt_median ~  Field + Sample.type * Season)
evals <- eigenvals(Ps_obj_ord1) # /sum(eigenvals(Ps_obj_ord)) * 100

Ps_obj_filt_median %>% 
  plot_ordination(., Ps_obj_ord1, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>% 
  mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>% 
  mutate_at(., "Sample.type", ~fct_relevel(., "Soil", "Green tea", "Rooibos")) %>% 
  dplyr::rename(., `Sample type` = Sample.type) ->
  ord_df

# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
#         width = 10, height = 7.2)

p_ord <- ggplot(ord_df,
             aes(
               x = CAP1,
               y = CAP2,
               shape = Field,
               color = `Sample type`
             )) +
  geom_point(size = 4, alpha = 2 / 3) +
  theme_bw(base_size = 14) +
  scale_colour_manual(values = pal_locuszoom("default")(3)[c(2, 3, 1)]) +
  stat_ellipse(
    aes(x = CAP1, y = CAP2, group = `Sample type`),
    alpha = 0.5,
    type = "norm",
    level = 0.95,
    linetype = 2,
    inherit.aes = FALSE
  ) +
  labs(x = sprintf("CAP1 [%s%%]", round(evals[1], 1)), 
       y = sprintf("CAP2 [%s%%]", round(evals[2], 1))) +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  facet_wrap(~ Season)

# # Now add the environmental variables as arrows
# arrowmat <- scores(Ps_obj_ord1, display = "bp") # bipplot arrows
# # Add labels, make a data.frame
# arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
# 
# arrowdf %<>%
#   mutate(labels = fct_recode(labels,
#     "Braunerde" = "FieldBraunerde",
#     "Kolluvisol" = "FieldKolluvisol",
#     "Rooibos" = "Sample.typeRooibos",
#     "Soil" = "Sample.typeSoil"
#   ))
# arrowdf %<>% dplyr::slice(., c(1:4))
# 
# # Define the arrow aesthetic mapping
# arrow_map = aes(xend = CAP1, yend = CAP2, x = 0, y = 0, shape = NULL, color = NULL, 
#     label = labels)
# label_map = aes(x = 1.2 * CAP1, y = 1.2 * CAP2, shape = NULL, color = NULL, 
#     label = labels)
# # Make a new graphic
# arrowhead = arrow(length = unit(0.05, "npc"))
# p_ord <- p_ord +
#   geom_segment(
#     arrow_map,
#     size = 0.5,
#     data = arrowdf,
#     color = "gray",
#     arrow = arrowhead
#   ) + 
#   geom_text(label_map, size = 2, data = arrowdf)
print(p_ord)

# dev.off()

# ggsave(
#   paste0(p_ord.file, ".png"),
#   p_ord,
#   device = "png",
#   width = 10,
#   height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))

Variance partitioning models and ordinations - soils excluded

Partitioning the data using discrete distance

Ps_obj_filt_median_tea <- subset_samples(Ps_obj_filt_median, Sample.type != "Soil")
adonis(
  otu_table(Ps_obj_filt_median_tea) ~ Lib.size,
  data = get_variable(Ps_obj_filt_median_tea),
  method = "horn",
  permutations = 9999
)

## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_tea) ~ Lib.size,      data = get_variable(Ps_obj_filt_median_tea), permutations = 9999,      method = "horn") 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Lib.size   1     0.292 0.29248 0.78766 0.00918  0.575
## Residuals 85    31.563 0.37133         0.99082       
## Total     86    31.855                 1.00000

adonis(
  otu_table(Ps_obj_filt_median_tea) ~ Field + Sample.type * Season,
  data = get_variable(Ps_obj_filt_median_tea),
  method = "horn",
  permutations = 999
)

## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_tea) ~ Field +      Sample.type * Season, data = get_variable(Ps_obj_filt_median_tea),      permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                    Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field               2     2.254  1.1270  7.4101 0.07076  0.001 ***
## Sample.type         1     2.608  2.6084 17.1497 0.08188  0.001 ***
## Season              3    11.925  3.9751 26.1360 0.37436  0.001 ***
## Sample.type:Season  3     3.356  1.1188  7.3558 0.10536  0.001 ***
## Residuals          77    11.711  0.1521         0.36764           
## Total              86    31.855                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# sample_data(Ps_obj_filt_median)$Type.season <- paste(sample_data(Ps_obj_filt)$Sample.type, sample_data(Ps_obj_filt)$Season) # too many comparisons

# Ps_obj_filt_median_s <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median) # for testing

# Compare sample-type pairs
mod_pairwise3 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median_tea),
  sample_data(Ps_obj_filt_median_tea)$Sample.type,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise3)

##                  pairs total.DF  F.Model         R2 p.value p.adjusted sig
## 1 Green tea vs Rooibos       86 7.952102 0.08555053   0.001      0.001  **

(sig_pairs3 <- as.character(mod_pairwise3$pairs[mod_pairwise3$p.adjusted < 0.05]))

## [1] "Green tea vs Rooibos"

# (meaningful_sig_pairs3 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))

# compare seasons
mod_pairwise4 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median_tea),
  sample_data(Ps_obj_filt_median_tea)$Season,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise4)

##              pairs total.DF   F.Model        R2 p.value p.adjusted sig
## 1 Winter vs Spring       43  8.533529 0.1688687   0.001     0.0012   *
## 2 Winter vs Summer       43 33.801961 0.4459246   0.001     0.0012   *
## 3 Winter vs Autumn       46 24.599278 0.3534416   0.001     0.0012   *
## 4 Spring vs Summer       39 20.786718 0.3535955   0.001     0.0012   *
## 5 Spring vs Autumn       42 13.464088 0.2472104   0.001     0.0012   *
## 6 Summer vs Autumn       42  3.578983 0.0802841   0.004     0.0040   *

(sig_pairs4 <- as.character(mod_pairwise4$pairs[mod_pairwise4$p.adjusted < 0.05]))

## [1] "Winter vs Spring" "Winter vs Summer" "Winter vs Autumn" "Spring vs Summer"
## [5] "Spring vs Autumn" "Summer vs Autumn"

# (meaningful_sig_pairs2 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))

Calculate ordinations

Ps_obj_ord3 <- ordinate(Ps_obj_filt_median_tea, "CAP", "horn", formula = Ps_obj_filt_median_tea ~  Field + Sample.type * Season)
evals <- eigenvals(Ps_obj_ord3) # /sum(eigenvals(Ps_obj_ord)) * 100

Ps_obj_filt_median_tea %>% 
  plot_ordination(., Ps_obj_ord3, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>% 
  mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>% 
  mutate_at(., "Sample.type", ~fct_relevel(., "Green tea", "Rooibos")) %>% 
  dplyr::rename(., `Sample type` = Sample.type) ->
  ord_df3

# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
#         width = 10, height = 7.2)

p_ord3 <- ggplot(ord_df3,
             aes(
               x = CAP1,
               y = CAP2,
               shape = Field,
               color = `Sample type`
             )) +
  geom_point(size = 4, alpha = 2 / 3) +
  theme_bw(base_size = 14) +
  scale_colour_manual(values = pal_locuszoom("default")(3)[c(3, 1)]) +
  stat_ellipse(
    aes(x = CAP1, y = CAP2, group = `Sample type`),
    alpha = 0.5,
    type = "norm",
    level = 0.95,
    linetype = 2,
    inherit.aes = FALSE
  ) +
  labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)), 
       y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  facet_wrap(~ Season)

print(p_ord3)

# dev.off()

# ggsave(
#   paste0(p_ord.file, ".png"),
#   p_ord,
#   device = "png",
#   width = 10,
#   height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))

Test differences between samples on the phylum and order levels

Taxa_tests_phylum1 <- STAMPR(Ps_obj_filt_median, 
                             vars2test = "Sample.type", 
                             sig_pairs = sig_pairs1, 
                             outputfile = paste0(Proj_name, "_Sample.type"))

plotSTAMPR(Taxa_tests_phylum1, pair = "Green tea vs Rooibos")

plotSTAMPR(Taxa_tests_phylum1, pair = "Green tea vs Soil")

plotSTAMPR(Taxa_tests_phylum1, pair = "Rooibos vs Soil")

Taxa_tests_order1 <- STAMPR(Ps_obj_filt_median, vars2test = "Sample.type", rank = "Order", sig_pairs = sig_pairs1, outputfile = paste0(Proj_name, "_Sample.type"))
plotSTAMPR(Taxa_tests_order1, pair = "Green tea vs Rooibos", tax_level = "Order")

plotSTAMPR(Taxa_tests_order1, pair = "Green tea vs Soil", tax_level = "Order")

plotSTAMPR(Taxa_tests_order1, pair = "Rooibos vs Soil", tax_level = "Order")

Taxa_tests_phylum2 <- STAMPR(Ps_obj_filt_median, vars2test = "Season", sig_pairs = sig_pairs2, outputfile = paste0(Proj_name, "_Season"))
plotSTAMPR(Taxa_tests_phylum2, pair = "Winter vs Summer")

Taxa_tests_order2 <- STAMPR(Ps_obj_filt_median, vars2test = "Season", rank = "Order", sig_pairs = sig_pairs2, outputfile = paste0(Proj_name, "_Season"))
plotSTAMPR(Taxa_tests_order2, pair = "Winter vs Summer", tax_level = "Order")

Taxa_tests_phylum4 <- STAMPR(Ps_obj_filt_median_tea, vars2test = "Season", sig_pairs = sig_pairs4, outputfile = paste0(Proj_name, "Teabags_Season"))
plotSTAMPR(Taxa_tests_phylum4, pair = "Winter vs Summer")

Taxa_tests_order4 <- STAMPR(Ps_obj_filt_median_tea, vars2test = "Season", rank = "Order", sig_pairs = sig_pairs4, outputfile = paste0(Proj_name, "Teabags_Season"))
plotSTAMPR(Taxa_tests_order4, pair = "Winter vs Summer", tax_level = "Order")

Differential abundance models

Tag rare phyla (for plotting purposes only)

Ps_obj_filt_media_glom <- tax_glom(Ps_obj_filt_median, 
                             "Order", 
                             NArm = TRUE) # glomerate to the Order level
Ps_obj_filt_media_glom_rel <- transform_sample_counts(Ps_obj_filt_media_glom, function(x) x / sum(x)) # transform to rel. ab.
Ps_obj_filt_media_glom_rel_DF <- psmelt(Ps_obj_filt_media_glom_rel) # generate a df
Ps_obj_filt_media_glom_rel_DF$Order %<>% as.character() # factor to char

# group dataframe by Order, calculate median rel. abundance
Ps_obj_filt_media_glom_rel_DF %>%
  group_by(Order) %>%
  summarise(median = median(Abundance)) ->
  medians

# find Phyla whose median rel. abund. is less than 0.5%
Rare_phyla <- medians[medians$median <= 0.005, ]$Order

# change their name to "Rare"
Ps_obj_filt_media_glom_rel_DF[Ps_obj_filt_media_glom_rel_DF$Order %in% Rare_phyla, ]$Order <- 'Rare'

# re-group
Ps_obj_filt_media_glom_rel_DF %>%
  group_by(Order) %>%
  summarise(Abundance = sum(Abundance)) %>% 
  arrange(desc(Abundance)) -> Taxa_rank

Detect differentially abundant OTUs using ALDEx2 (Fernandes et al. 2013)

Sample type differences

significance = 0.05

# run full model 
data2test <- t(otu_table(Ps_obj_filt_median))
# comparison <- as.character(get_variable(Ps_obj_filt_median, "Sample.type"))
#ALDEx_full <- aldex.clr(data2test, comparison, mc.samples = 128, denom = "iqlr", verbose = TRUE, useMC = TRUE) # iqlr for slight assymetry in composition
#ALDEx_full_glm <- aldex.glm(ALDEx_full, comparison, useMC = TRUE) # for more than two conditions
#sig_taxa <- rownames(ALDEx_full_glm)[ALDEx_full_glm$glm.eBH < 0.05] # save names of taxa that are significant under the full model
#write.csv(sig_taxa, file = "Aldex_full_significant_taxa.csv")

# Pairwise comparisons
# 
ALDEx_comparisons <- list()
ALDEx_comparisons$Comparisons <- sig_pairs1

# Ps_obj_filt_median_subset <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median)

for (j in seq(1, length(ALDEx_comparisons$Comparisons))) {
  # print(j)
  ALDEx_comparisons$Comparisons[j] %>% 
    str_split(., " vs ", simplify = FALSE) %>% 
    unlist() ->
    comparison_string
  Ps_obj_filt_median %>%
    subset_samples(Sample.type == comparison_string[1] | Sample.type == comparison_string[2]) ->
    # subset_samples(Treatment == comparison_string[2] | Treatment == comparison_string[4]) %>%
    # subset_samples(Year == "2016" | Year == "2017" | Year == "2018") ->
    Ps_obj_filt_median_pairwise

#  Remove species with prevalence < 10%
  Ps_obj_filt_median_pairwise_s <- DropRareSpecies(Ps_obj = Ps_obj_filt_median_pairwise, prevalence = 0.1)
    
  sample_data(Ps_obj_filt_median_pairwise_s)$Sample.type %<>% fct_relevel(., str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1])
  
  # make Joint.sample.name for matching OTUs between compared samples (for GGPlotTopOTUs)
  suppressWarnings(
  sample_data(Ps_obj_filt_median_pairwise_s) %<>%
    as(., "data.frame") %>% 
    rownames_to_column() %>% 
    mutate_if(is.factor, as.character) %>% 
    mutate(Joint.sample.name  = paste0(.$Sample.type, "_", .$Field, "_", .$Season, "_", .$Replicate)) %>% 
    mutate_at("Joint.sample.name", ~str_replace_all(., paste0(unique(comparison_string), collapse = "|"), "")) %>% # remove the levels participating in the comparison from the names
    mutate_at("Joint.sample.name", ~str_replace_all(., "^_", "")) %>% # remove first "_"
    column_to_rownames()
  )
    
  ALDEx2plot_pairwise <- CalcALDEx(
    physeq_obj = Ps_obj_filt_median_pairwise_s,
    vars2test = "Sample.type",
    rare_phyla = Rare_phyla,
    sig_level = significance,
    LFC = 0
    )
  ALDEx_comparisons$Results[[j]] <- ALDEx2plot_pairwise # store results
  ALDEx_comparisons$Results[[j]]$Var1 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1]
  ALDEx_comparisons$Results[[j]]$Var2 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[2]
  ALDEX_summary <- tibble(Label = c(paste0("⬆", sum( ALDEx_comparisons$Results[[j]]$effect > 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " ⬇", sum( ALDEx_comparisons$Results[[j]]$effect < 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " (", nrow( ALDEx_comparisons$Results[[j]]), ")")))
  
  ALDEx2plot_pairwise %>%
    filter(Significance == "Pass") %>%
    dplyr::select(OTU, baseMean, effect, Order, Class, Order, Family, Genus) %>%
    arrange(desc(abs(effect))) ->
    ALDEx2plot_pairwise_results
    
  # print(ALDEx2plot_pairwise_results %>% 
  #   kable(., digits = c(2), caption = "Significantly different taxa:") %>%
  #   kable_styling(
  #     bootstrap_options = c("striped", "hover", "condensed", "responsive"),
  #     full_width = F
  #   ))
  
  write.csv(ALDEx2plot_pairwise_results, file = paste0("Aldex", "_", paste0(comparison_string, collapse = "_"), ".csv"))
  
  # Plot ALDEX plot
  p1 <- GGPlotALDExTax(ALDEx2plot_pairwise, OTU_labels = FALSE, Taxa = "Order", sig_level = significance) +
    # ggtitle(ALDEx_comparisons$Comparisons[j]) +
    geom_text(
    data    = ALDEX_summary,
    mapping = aes(x = Inf, y = Inf, label = Label),
    hjust   = 1.1,
    vjust   = 1.6
  ) 
  p1 <- p1 + labs(title = ALDEx_comparisons$Comparisons[j])
  print(p1)
  
  # Plot OTU plots
  # GGPlotOTU(Ps_obj_filt_subset_pairwise_s, vars2test = "Spill.Treatment", "Seq_8")
  p2 <- GGPlotTopOTUs(
    physeq_obj = Ps_obj_filt_median_pairwise_s,
    vars2test = "Sample.type",
    ALDEx_obj = ALDEx2plot_pairwise,
    tax_level = "Order",
    rank_by = "effect",
    Ntop = 12
  )
  print(p2)
}

phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 87 samples ] sample_data() Sample Data: [ 87 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 266 taxa and 87 samples ] sample_data() Sample Data: [ 87 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 266 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 75 samples ] sample_data() Sample Data: [ 75 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 275 taxa and 75 samples ] sample_data() Sample Data: [ 75 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 275 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 72 samples ] sample_data() Sample Data: [ 72 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 297 taxa and 72 samples ] sample_data() Sample Data: [ 72 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 297 taxa by 6 taxonomic ranks ]

filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU[common_taxa <- filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU %in% filter(ALDEx_comparisons$Results[[2]], Significance == "Pass")$OTU]

[1] "Seq_41" "Seq_42" "Seq_78" "Seq_204"

Seasonal differences

Comparing the seasonsal differences in teabag samples only

significance = 0.05

data2test <- t(otu_table(Ps_obj_filt_median_tea))

# Pairwise comparisons
ALDEx_comparisons <- list()
ALDEx_comparisons$Comparisons <- sig_pairs2

# Ps_obj_filt_median_tea_subset <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median_tea), TRUE)[1:100]), Ps_obj_filt_median_tea)

for (j in seq(1, length(ALDEx_comparisons$Comparisons))) {
  # print(j)
  ALDEx_comparisons$Comparisons[j] %>% 
    str_split(., " vs ", simplify = FALSE) %>% 
    unlist() ->
    comparison_string
  Ps_obj_filt_median_tea %>%
    subset_samples(Season == comparison_string[1] | Season == comparison_string[2]) ->
    # subset_samples(Treatment == comparison_string[2] | Treatment == comparison_string[4]) %>%
    # subset_samples(Year == "2016" | Year == "2017" | Year == "2018") ->
    Ps_obj_filt_median_tea_pairwise

#  Remove species with prevalence < 10%
  Ps_obj_filt_median_tea_pairwise_s <- DropRareSpecies(Ps_obj = Ps_obj_filt_median_tea_pairwise, prevalence = 0.1)
    
  sample_data(Ps_obj_filt_median_tea_pairwise_s)$Season %<>% fct_relevel(., str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1])
  
  # make Joint.sample.name for matching OTUs between compared samples (for GGPlotTopOTUs)
  suppressWarnings(
  sample_data(Ps_obj_filt_median_tea_pairwise_s) %<>% 
    as(., "data.frame") %>% 
    rownames_to_column() %>% 
    mutate_if(is.factor, as.character) %>% 
    mutate(Joint.sample.name  = paste0(.$Sample.type, "_", .$Field, "_", .$Season, "_", .$Replicate)) %>% 
    mutate_at("Joint.sample.name", ~str_replace_all(., paste0(unique(comparison_string), collapse = "|"), "")) %>% # remove the levels participating in the comparison from the names
    mutate_at("Joint.sample.name", ~str_replace_all(., "^_", "")) %>% # remove first "_"
    mutate_at("Joint.sample.name", ~str_replace_all(., "__", "_")) %>% # remove double "_"
    column_to_rownames()
  )
    
  ALDEx2plot_pairwise <- CalcALDEx(
    physeq_obj = Ps_obj_filt_median_tea_pairwise_s,
    vars2test = "Season",
    rare_phyla = Rare_phyla,
    sig_level = significance,
    LFC = 0
    )
  ALDEx_comparisons$Results[[j]] <- ALDEx2plot_pairwise # store results
  ALDEx_comparisons$Results[[j]]$Var1 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1]
  ALDEx_comparisons$Results[[j]]$Var2 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[2]
  ALDEX_summary <- tibble(Label = c(paste0("⬆", sum( ALDEx_comparisons$Results[[j]]$effect > 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " ⬇", sum( ALDEx_comparisons$Results[[j]]$effect < 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " (", nrow( ALDEx_comparisons$Results[[j]]), ")")))
  
  ALDEx2plot_pairwise %>%
    filter(Significance == "Pass") %>%
    dplyr::select(OTU, baseMean, effect, Order, Class, Order, Family, Genus) %>%
    arrange(desc(abs(effect))) ->
    ALDEx2plot_pairwise_results
    
  # print(ALDEx2plot_pairwise_results %>% 
  #   kable(., digits = c(2), caption = "Significantly different taxa:") %>%
  #   kable_styling(
  #     bootstrap_options = c("striped", "hover", "condensed", "responsive"),
  #     full_width = F
  #   ))
  
  write.csv(ALDEx2plot_pairwise_results, file = paste0("Aldex", "_", paste0(comparison_string, collapse = "_"), ".csv"))
  
  # Plot ALDEX plot
  p1 <- GGPlotALDExTax(ALDEx2plot_pairwise, OTU_labels = FALSE, sig_level = significance) +
    # ggtitle(ALDEx_comparisons$Comparisons[j]) +
    geom_text(
    data    = ALDEX_summary,
    mapping = aes(x = Inf, y = Inf, label = Label),
    hjust   = 1.1,
    vjust   = 1.6
  ) 
  p1 <- p1 + labs(title = ALDEx_comparisons$Comparisons[j])
  print(p1)
  
  # Plot OTU plots
  p2 <- GGPlotTopOTUs(
    physeq_obj = Ps_obj_filt_median_tea_pairwise_s,
    vars2test = "Season",
    ALDEx_obj = ALDEx2plot_pairwise,
    tax_level = "Order",
    rank_by = "effect",
    Ntop = 12
  )
  print(p2)
}

phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 262 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 262 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 293 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 293 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 316 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 316 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 40 samples ] sample_data() Sample Data: [ 40 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 196 taxa and 40 samples ] sample_data() Sample Data: [ 40 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 196 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 43 samples ] sample_data() Sample Data: [ 43 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 198 taxa and 43 samples ] sample_data() Sample Data: [ 43 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 198 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 402 taxa and 43 samples ] sample_data() Sample Data: [ 43 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ] phyloseq-class experiment-level object otu_table() OTU Table: [ 186 taxa and 43 samples ] sample_data() Sample Data: [ 43 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 186 taxa by 6 taxonomic ranks ]

filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU[common_taxa <- filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU %in% filter(ALDEx_comparisons$Results[[2]], Significance == "Pass")$OTU]

[1] "Seq_3" "Seq_17" "Seq_20" "Seq_27" "Seq_31" "Seq_33" "Seq_36" "Seq_40" [9] "Seq_43" "Seq_44" "Seq_50" "Seq_55" "Seq_57" "Seq_64" "Seq_67" "Seq_70" [17] "Seq_72" "Seq_74" "Seq_80" "Seq_85" "Seq_95" "Seq_107" "Seq_128" "Seq_203" [25] "Seq_240"

sessioninfo::session_info() %>%
  details::details(
    summary = 'Current session info',
    open    = TRUE
 )

Current session info

─ Session info ─────────────────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 4.0.3 (2020-10-10)
 os       Ubuntu 18.04.5 LTS          
 system   x86_64, linux-gnu           
 ui       X11                         
 language en_GB                       
 collate  en_DK.UTF-8                 
 ctype    en_DK.UTF-8                 
 tz       Europe/Prague               
 date     2021-05-21                  

─ Packages ─────────────────────────────────────────────────────────────────────────────
 package              * version  date       lib source                                  
 ade4                 * 1.7-16   2020-10-28 [1] CRAN (R 4.0.2)                          
 affy                   1.68.0   2020-10-27 [1] Bioconductor                            
 affyio                 1.60.0   2020-10-27 [1] Bioconductor                            
 ALDEx2               * 1.22.0   2020-10-27 [1] Bioconductor                            
 ape                    5.5      2021-04-25 [1] CRAN (R 4.0.3)                          
 assertthat             0.2.1    2019-03-21 [1] CRAN (R 4.0.2)                          
 backports              1.2.1    2020-12-09 [1] CRAN (R 4.0.2)                          
 Biobase              * 2.50.0   2020-10-27 [1] Bioconductor                            
 BiocGenerics         * 0.36.1   2021-04-16 [1] Bioconductor                            
 BiocManager            1.30.15  2021-05-11 [1] CRAN (R 4.0.3)                          
 BiocParallel           1.24.1   2020-11-06 [1] Bioconductor                            
 biomformat             1.18.0   2020-10-27 [1] Bioconductor                            
 Biostrings             2.58.0   2020-10-27 [1] Bioconductor                            
 bitops                 1.0-7    2021-04-24 [1] CRAN (R 4.0.3)                          
 broom                * 0.7.6    2021-04-05 [1] CRAN (R 4.0.3)                          
 bslib                  0.2.5.1  2021-05-18 [1] CRAN (R 4.0.3)                          
 cellranger             1.1.0    2016-07-27 [1] CRAN (R 4.0.2)                          
 cli                    2.5.0    2021-04-26 [1] CRAN (R 4.0.3)                          
 clipr                  0.7.1    2020-10-08 [1] CRAN (R 4.0.2)                          
 cluster                2.1.2    2021-04-17 [1] CRAN (R 4.0.3)                          
 codetools              0.2-18   2020-11-04 [1] CRAN (R 4.0.2)                          
 colorspace             2.0-1    2021-05-04 [1] CRAN (R 4.0.3)                          
 cowplot              * 1.1.1    2020-12-30 [1] CRAN (R 4.0.2)                          
 crayon                 1.4.1    2021-02-08 [1] CRAN (R 4.0.3)                          
 data.table             1.14.0   2021-02-21 [1] CRAN (R 4.0.3)                          
 DBI                    1.1.1    2021-01-15 [1] CRAN (R 4.0.3)                          
 dbplyr                 2.1.1    2021-04-06 [1] CRAN (R 4.0.3)                          
 DelayedArray           0.16.3   2021-03-24 [1] Bioconductor                            
 desc                   1.3.0    2021-03-05 [1] CRAN (R 4.0.3)                          
 details                0.2.1    2020-01-12 [1] CRAN (R 4.0.2)                          
 digest                 0.6.27   2020-10-24 [1] CRAN (R 4.0.2)                          
 doParallel           * 1.0.16   2020-10-16 [1] CRAN (R 4.0.2)                          
 dplyr                * 1.0.6    2021-05-05 [1] CRAN (R 4.0.3)                          
 ellipsis               0.3.2    2021-04-29 [1] CRAN (R 4.0.3)                          
 evaluate               0.14     2019-05-28 [1] CRAN (R 4.0.2)                          
 extrafont            * 0.17     2014-12-08 [1] CRAN (R 4.0.2)                          
 extrafontdb            1.0      2012-06-11 [1] CRAN (R 4.0.2)                          
 fansi                  0.4.2    2021-01-15 [1] CRAN (R 4.0.3)                          
 farver                 2.1.0    2021-02-28 [1] CRAN (R 4.0.3)                          
 forcats              * 0.5.1    2021-01-27 [1] CRAN (R 4.0.3)                          
 foreach              * 1.5.1    2020-10-15 [1] CRAN (R 4.0.2)                          
 fs                     1.5.0    2020-07-31 [1] CRAN (R 4.0.2)                          
 generics               0.1.0    2020-10-31 [1] CRAN (R 4.0.2)                          
 GenomeInfoDb           1.26.7   2021-04-08 [1] Bioconductor                            
 GenomeInfoDbData       1.2.4    2021-05-05 [1] Bioconductor                            
 GenomicRanges          1.42.0   2020-10-27 [1] Bioconductor                            
 ggplot2              * 3.3.3    2020-12-30 [1] CRAN (R 4.0.2)                          
 ggpomological        * 0.1.2    2020-08-13 [1] Github (gadenbuie/ggpomological@69f3815)
 ggrepel              * 0.9.1    2021-01-15 [1] CRAN (R 4.0.3)                          
 ggsci                * 2.9      2018-05-14 [1] CRAN (R 4.0.2)                          
 ggthemes             * 4.2.4    2021-01-20 [1] CRAN (R 4.0.3)                          
 glue                   1.4.2    2020-08-27 [1] CRAN (R 4.0.2)                          
 gridExtra            * 2.3      2017-09-09 [1] CRAN (R 4.0.2)                          
 gtable                 0.3.0    2019-03-25 [1] CRAN (R 4.0.2)                          
 haven                  2.4.1    2021-04-23 [1] CRAN (R 4.0.3)                          
 hexbin                 1.28.2   2021-01-08 [1] CRAN (R 4.0.2)                          
 highr                  0.9      2021-04-16 [1] CRAN (R 4.0.3)                          
 hms                    1.1.0    2021-05-17 [1] CRAN (R 4.0.3)                          
 htmltools              0.5.1.1  2021-01-22 [1] CRAN (R 4.0.3)                          
 httr                   1.4.2    2020-07-20 [1] CRAN (R 4.0.2)                          
 igraph                 1.2.6    2020-10-06 [1] CRAN (R 4.0.2)                          
 IRanges                2.24.1   2020-12-12 [1] Bioconductor                            
 iterators            * 1.0.13   2020-10-15 [1] CRAN (R 4.0.2)                          
 jquerylib              0.1.4    2021-04-26 [1] CRAN (R 4.0.3)                          
 jsonlite               1.7.2    2020-12-09 [1] CRAN (R 4.0.2)                          
 kableExtra           * 1.3.4    2021-02-20 [1] CRAN (R 4.0.3)                          
 knitr                * 1.33     2021-04-24 [1] CRAN (R 4.0.3)                          
 labeling               0.4.2    2020-10-20 [1] CRAN (R 4.0.2)                          
 lattice              * 0.20-44  2021-05-02 [1] CRAN (R 4.0.3)                          
 lifecycle              1.0.0    2021-02-15 [1] CRAN (R 4.0.3)                          
 limma                  3.46.0   2020-10-27 [1] Bioconductor                            
 lubridate              1.7.10   2021-02-26 [1] CRAN (R 4.0.3)                          
 magrittr             * 2.0.1    2020-11-17 [1] CRAN (R 4.0.2)                          
 MASS                 * 7.3-54   2021-05-03 [1] CRAN (R 4.0.3)                          
 Matrix                 1.3-3    2021-05-04 [1] CRAN (R 4.0.3)                          
 MatrixGenerics         1.2.1    2021-01-30 [1] Bioconductor                            
 matrixStats            0.58.0   2021-01-29 [1] CRAN (R 4.0.3)                          
 mgcv                   1.8-35   2021-04-18 [1] CRAN (R 4.0.3)                          
 modelr                 0.1.8    2020-05-19 [1] CRAN (R 4.0.2)                          
 multtest               2.46.0   2020-10-27 [1] Bioconductor                            
 munsell                0.5.0    2018-06-12 [1] CRAN (R 4.0.2)                          
 NADA                 * 1.6-1.1  2020-03-22 [1] CRAN (R 4.0.2)                          
 nlme                   3.1-152  2021-02-04 [1] CRAN (R 4.0.3)                          
 permute              * 0.9-5    2019-03-12 [1] CRAN (R 4.0.2)                          
 phyloseq             * 1.34.0   2020-10-27 [1] Bioconductor                            
 pillar                 1.6.1    2021-05-16 [1] CRAN (R 4.0.3)                          
 pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 4.0.2)                          
 plyr                   1.8.6    2020-03-03 [1] CRAN (R 4.0.2)                          
 png                    0.1-7    2013-12-03 [1] CRAN (R 4.0.2)                          
 preprocessCore         1.52.1   2021-01-08 [1] Bioconductor                            
 prettyunits            1.1.1    2020-01-24 [1] CRAN (R 4.0.2)                          
 progress               1.2.2    2019-05-16 [1] CRAN (R 4.0.2)                          
 ps                     1.6.0    2021-02-28 [1] CRAN (R 4.0.3)                          
 purrr                * 0.3.4    2020-04-17 [1] CRAN (R 4.0.2)                          
 R6                     2.5.0    2020-10-28 [1] CRAN (R 4.0.2)                          
 ragg                 * 1.1.2    2021-03-17 [1] CRAN (R 4.0.3)                          
 Rcpp                   1.0.6    2021-01-15 [1] CRAN (R 4.0.3)                          
 RCurl                  1.98-1.3 2021-03-16 [1] CRAN (R 4.0.3)                          
 readr                * 1.4.0    2020-10-05 [1] CRAN (R 4.0.2)                          
 readxl                 1.3.1    2019-03-13 [1] CRAN (R 4.0.2)                          
 reprex                 2.0.0    2021-04-02 [1] CRAN (R 4.0.3)                          
 reshape2               1.4.4    2020-04-09 [1] CRAN (R 4.0.2)                          
 rhdf5                  2.34.0   2020-10-27 [1] Bioconductor                            
 rhdf5filters           1.2.1    2021-05-03 [1] Bioconductor                            
 Rhdf5lib               1.12.1   2021-01-26 [1] Bioconductor                            
 rlang                  0.4.11   2021-04-30 [1] CRAN (R 4.0.3)                          
 rmarkdown            * 2.8      2021-05-07 [1] CRAN (R 4.0.3)                          
 rprojroot              2.0.2    2020-11-15 [1] CRAN (R 4.0.2)                          
 rstudioapi             0.13     2020-11-12 [1] CRAN (R 4.0.2)                          
 Rttf2pt1               1.3.8    2020-01-10 [1] CRAN (R 4.0.2)                          
 rvest                  1.0.0    2021-03-09 [1] CRAN (R 4.0.3)                          
 S4Vectors              0.28.1   2020-12-09 [1] Bioconductor                            
 sass                   0.4.0    2021-05-12 [1] CRAN (R 4.0.3)                          
 scales               * 1.1.1    2020-05-11 [1] CRAN (R 4.0.2)                          
 sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 4.0.2)                          
 stringi                1.6.2    2021-05-17 [1] CRAN (R 4.0.3)                          
 stringr              * 1.4.0    2019-02-10 [1] CRAN (R 4.0.2)                          
 SummarizedExperiment   1.20.0   2020-10-27 [1] Bioconductor                            
 survival             * 3.2-11   2021-04-26 [1] CRAN (R 4.0.3)                          
 svglite              * 2.0.0    2021-02-20 [1] CRAN (R 4.0.3)                          
 systemfonts            1.0.2    2021-05-11 [1] CRAN (R 4.0.3)                          
 textshaping            0.3.4    2021-05-11 [1] CRAN (R 4.0.3)                          
 tibble               * 3.1.2    2021-05-16 [1] CRAN (R 4.0.3)                          
 tidyr                * 1.1.3    2021-03-03 [1] CRAN (R 4.0.3)                          
 tidyselect             1.1.1    2021-04-30 [1] CRAN (R 4.0.3)                          
 tidyverse            * 1.3.1    2021-04-15 [1] CRAN (R 4.0.3)                          
 truncnorm            * 1.0-8    2018-02-27 [1] CRAN (R 4.0.2)                          
 utf8                   1.2.1    2021-03-12 [1] CRAN (R 4.0.3)                          
 vctrs                  0.3.8    2021-04-29 [1] CRAN (R 4.0.3)                          
 vegan                * 2.5-7    2020-11-28 [1] CRAN (R 4.0.3)                          
 viridisLite            0.4.0    2021-04-13 [1] CRAN (R 4.0.3)                          
 vsn                  * 3.58.0   2020-10-27 [1] Bioconductor                            
 webshot                0.5.2    2019-11-22 [1] CRAN (R 4.0.2)                          
 withr                  2.4.2    2021-04-18 [1] CRAN (R 4.0.3)                          
 xfun                   0.23     2021-05-15 [1] CRAN (R 4.0.3)                          
 xml2                   1.3.2    2020-04-23 [1] CRAN (R 4.0.2)                          
 XVector                0.30.0   2020-10-27 [1] Bioconductor                            
 yaml                   2.2.1    2020-02-01 [1] CRAN (R 4.0.2)                          
 zCompositions        * 1.3.4    2020-03-04 [1] CRAN (R 4.0.2)                          
 zlibbioc               1.36.0   2020-10-27 [1] Bioconductor                            

[1] /home/angel/R/library
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library

References

Fernandes AD, Macklaim JM, Linn TG et al. ANOVA-Like Differential Expression (ALDEx) Analysis for Mixed Population RNA-Seq. PLOS ONE 2013;8:e67019.